#include<iostream>
using namespace std;
#include<ctime>
#include<omp.h>

/*
 * Simpson：辛普生积分
 * 给定积分区间[a,b]，区间分割为n份，被积函数为func，返回积分值
 */
double Simpson(double a, double b, long long n, double (*func)(double)){

        double y1,y2,y,h;
        long long i;

        h = (b-a)/n;//步长

        y1 = 0;
	#pragma omp parallel for reduction(+:y1)
        for(i=0;i<n;i++){
                y1 += func(a+ (i+0.5)*h);
        }
        y2 = 0;
	#pragma omp parallel for reduction(+:y2)
        for(i=1;i<n;i++){
                y2 += func(a+ i*h);
        }
        y = 4*y1 + 2*y2 + func(a) + func(b);
        y *= h/6;
        return y;
}

double f(double x){
	return x*x;
}

int main(){

	double a=0,b=1;
	long long n = 1E9;

	time_t t_start = time(NULL);
	cout<<"integ = "<<Simpson(a,b,n,f)<<endl;
	time_t t_end = time(NULL);
	cout<<"Elapsed time: "<<(t_end-t_start)<<" s.\n";

	return 0;
}

/*
 * n = 1E9	serial = 12 s	parallel(6 core) = 2 s
 */
